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Abstract. Inferring network topology from dynamical observations is a fundamental 
problem pervading research on complex systems. Here, we present a simple, direct 
^ method to infer the structural connection topology of a network, given an observation 

of one collective dynamical trajectory. The general theoretical framework is applicable 
to arbitrary network dynamical systems described by ordinary differential equations. 
I— I No interference (external driving) is required and the type of dynamics is not 

^J restricted in any way. In particular, the observed dynamics may be arbitrarily 

Q^ complex; stationary, invariant or transient; synchronous or asynchronous and chaotic 

or periodic. Presupposing a knowledge of the functional form of the dynamical units 
and of the coupling functions between them, we present an analytical solution to the 
inverse problem of finding the network topology. Robust reconstruction is achieved 
in any sufficiently long generic observation of the system. We extend our method 
to simultaneously reconstruct both the entire network topology and all parameters 
^ appearing linear in the system's equations of motion. Reconstruction of network 

topology and system parameters is viable even in the presence of substantial external 



\^ noise. 

PACS numbers: 89.75.-k, 05.45.Xt 05.45.Tp 
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1. Background 

Understanding the relations between network topology and its collective dynamics is 
at the heart of current interdisciplinary research on networked systems p!]. Often it 
is possible to observe the dynamics of the individual units of the network, whereas 
the coupling strengths between them and the underlying network topology cannot be 
directly measured. Hence, various methods have been proposed to solve the inverse 
problem of inferring network structure from observation and control of dynamics. 

Perturbing a fixed point of a network dynamical system constitutes the simplest 
controlled intervention of a system. The method of Tegner et al. [2] perturbs the 
steady state expression levels of selected genes. Their iterative algorithm can reveal the 
structure of an underlying gene regulatory network by analysing resultant dynamical 
changes in the pattern of gene expression levels. A similar iterative method based on 
multiple regression, coupled with transcriptional perturbations to the fixed points of 
a genetic network, has been used to successfully identify a nine-gene sub-network [3]. 
A method introduced by Timme |4] extends reconstruction to networks of smoothly 
coupled limit-cycle oscillators with periodic collective dynamics. The underlying idea 
is that the asymptotic response dynamics of a network to different externally induced 
driving conditions is a function of its topology and of the (external) driving signals. Thus 
measuring the response to suitable driving signals in different experiments restricts the 
set of network topologies that are consistent with the driving-response pairs, yielding 
the network's topology for sufficiently many experiments. 

Can we reconstruct a network displaying collective dynamics richer than simple 
fixed points or limit-cycles? Yu et al. [5j introduced a synchronization method to identify 
networks of chaotic Lorenz oscillators up to N=17 units. Assuming full knowledge of all 
model parameters, the network topology of a clone of the system is varied progressively 
via error minimisation until it synchronizes with the original system. The topology 
of the clone is then recognised as that of the original network. An extension of 
this synchronisation method [6j involves additional 'control signals' to externally drive 
the system to steady states, allowing the inference of interaction topology for sparse, 
symmetric networks. 

Here, we introduce a simple, direct and intervention-free method to reconstruct 
networks of arbitrary topology from the mere observation of generic collective dynamics. 
Given the functional form of the intrinsic and interaction dynamics, we show that all 
other factors of the network dynamical system, such as network topology and coupling 
strengths (and even typical parameters), can be reliably and efficiently reconstructed. 

2. Theory of Direct Reconstruction from Dynamical Trajectories 

Given an observation of the collective trajectory of a dynamical network, how can we 
infer its underlying topology? We consider a dynamical system consisting of A^ units, 
where the dynamics of each unit is specified by an arbitrary set of dynamical equations. 
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Figure 1. Reconstruction of networks is possible for various types of collective 
dynamics. A 32-unit random network of Rossler oscillators ^ with fixed connection 
probability p = 0.5 exhibiting (a) a synchronous periodic state (a^ = 0.2, bi = 1.7 
and Ci = 4.0 for all units); (b) a synchronous chaotic state (a^ = 0.2, bi = 1.7 and 
Ci = 13.0 for all units); (c) an asynchronous periodic state (parameters randomly 
drawn independently from a^ G [0.10,0.101], bi e [0.10,0.101] and q G {4,6,12}); 
(d) an asynchronous chaotic state (parameters randomly drawn independently from 
ai e [0.0,0.38], bi e [0.0,2.0] and q G [13,16]) The x-values of three out of TV = 32 
oscillators are shown. Close to perfect reconstructions from all these distinct dynamics 
are shown in Fig. [2] 



and the interactions between the units take place via edges in the network. The units 
of the dynamical system are coupled on a directed graph of unknown connectivity with 
their dynamics satisfying 
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where i^j G {1, 2, . . . A^}, and Xi = 

i-th unit, and the functions f^, g^ : M^ -^ M^ mediate intrinsic and interaction dynamics 
of the D-dimensional unit, and are known. Methods based on copy-synchronisation ^]6\ 
rely on the construction of a new network, with dynamics governed by Eq.[T]and network 
parameters J^j that are tuned to that of the real network by an error minimisation 
procedure. Here, we reduce the same reconstruction problem by evaluating the states 
and their derivatives directly, recognising that the only remaining unknowns in Eq. [T] 
are the coupling strengths, which are to be determined. 

The dynamics of the d-th dimension of the i-th unit is given by 
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where t^ G M are the times we evaluate Equation [T] and we now write z for the rate of 
change ^z of a scalar variable z. If there are M such times, m G {1, . . . , M}, we have 



Inferring Network Topology from Complex Dynamics 4 

M equations of the form 

N 
Ad) _ Ad) , V- 7 Ad) (n^ 

for each dimension d of the local dynamical systems f- separately. As the equations 
([3]) are uncoupled for any two different dimensions d and d\ we treat these separately 
and drop the index (d) from now on. 

Repeated evaluations of the equations of motion ^ of the system at different times 
tm thus comprise a simple and implicit restriction on the network topology Jij as follows: 
writing X^^^ = Xi^m — fi.m-, these equations constitute the matrix equation 

X, = J,G, (4) 

where Xi G R^""^ , Ji G R^""^ and Q G R^""^ . Here the elements of the i-th row of J 
are given by Ji and comprise the sequence {Jij)je{i,...,N} of all input coupling strengths 
to unit i 

Can we rewrite this equation explicitly for Ji? Generically, M > N ^ and we wish 
to solve this overdetermined problem by minimising the error function given by 

M / N \ 2 

EiiJij = / ^ I Xijji — y ^ JikQikm I (5) 

m=l \ k=l ) 

for the best (in Euclidean (^2) norm) solution Ji. Here J^/e represents the reconstructed 
value of the real coupling strength Ji].. Equating to zero the derivatives of the error 
function with respect to the matrix elements, Qj—Ei I Jij =0, yields an analytical 
solution to £2 error-minimisation given by 

Ji = XGj{GiGjy' (6) 

and thus the set of input coupling strengths (and input connectivity) of unit i. 
Evaluating such equations for all i G {1,...,A^} yields the complete reconstructed 
network J. This mathematical form of minimum ^2-norm solution is implemented in 
many mathematical packages (e.g. as the mrdivide function in MATLAB [7]). 

3. Performance for Different Collective Dynamics 

How does this theoretical method perform in applications on data? To illustrate the 
performance of the method and its insensitivity to the type of dynamics, we apply it 
to four distinct collective dynamical states ranging from simple periodic synchronous 
dynamics to very complex, highly chaotic asynchronous states. We choose networks of 
Rossler oscillators that can exhibit a rich repertoire of collective dynamics from multi- 
dimensional chaos to periodicity and from global synchrony to asynchrony (see Fig. fl]), 
depending on local unit parameters and coupling functions. 
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Figure 2. Successful reconstruction of the topology of networks from very different 
dynamics. Each of the panels (a)-(d) shows a reconstruction of the network's topology 
from the dynamics shown in the respective panels (a)-(d) of Figure 1. The smaller 
panels show absolute differences (magnified by a factor of 10^) from the actual network 
(picked randomly from an ensemble of networks with connection probability p = 0.5). 
The reconstructions (a)-(d) used trajectories from very different dynamical states, as 
shown in Figure [l] (a)-(d). Reconstructions in (c) and (d) have lower errors as they 
utilise whole dynamical trajectories, instead of transients towards synchrony, as in (a) 
and (b). 



3.1. Successful reconstruction from different dynamics 

The dynamics of each Rossler oscillator [8j is given by the three ordinary differential 
equations 

N 



Jbi 






yi = Xi + aiyi 



(7) 



where a^, 6^, q are local parameters. The coupling function was set to / (x^, Xj) = Xj —Xi 
to induce synchronisation and to f{xi^Xj) = sin(xj) to prevent it. The local unit 
parameters a^, bi and q of the Rossler oscillators are chosen to induce either chaotic 
or periodic dynamics. The parameters were treated as unknown and are not needed 
for the reconstruction of the network topology. This is the case in this example, since 
the equations where the parameters of network topology occur (in the x-dimension) do 
not contain any other (unknown) parameter. We now demonstrate a reconstruction of 
the network in all four dynamical paradigms illustrated in Fig. [TJ periodic synchronous, 
chaotic synchronous, periodic asynchronous and chaotic asynchronous collective states. 
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Reconstruction in praxis works as follows: for networks exhibiting synchronised 
dynamics, the coupling function / (x^, Xj) = Xj — Xi is uniformly zero for all time for all 
units, revealing no information about the network topology. Nevertheless, the network 
can still be reconstructed from its transient dynamics towards the synchronous state. In 
general, by substituting X^,^ = Xi {tm)+yi (tm)+^i (^m) and Q^- ^ = / (^^(tm), ^j(^m)) in 
^ , we find the least-squares reconstructed network according to ^ . Since each unit has 
at most A^ — 1 incoming links of unknown weights, the diagonals in the reconstructed 
J are zero. A reliable reconstruction of the network from trajectories (as shown in 
Figure [l]) is illustrated in Figure [2| for all four dynamical paradigms considered. 

3.2. Quality of Reconstruction 

How accurate is such a reconstruction? The quality of reconstruction 

is defined as the fraction of coupling strengths that are considered correct. Here a <1 
is the required accuracy of the coupling strengths and H is the Heaviside step function, 
H [x) = 1 for X > and is H{x) = otherwise. The normalized element-wise difference 
between the reconstructed and real network is 



LXJij . Jij J ij 



I x^^max) \^ ) 



where Jmax = max^/j/ < \Ji'j'\ , Ji'jn \- Typically, the quality of reconstruction increases 
with both the sampling rate and the length of trajectory observed, becoming close to 1 
even at lower sampling rates for longer times (see Fig. [s]). 

3.3. Required observation time 

What is the minimum length of observation required to reconstruct a network? We 
define 

T,,,,^:=min{r|g,(r)>g} (10) 

to be the minimum length of time of observation at a sampling rate u required for 
accurate reconstruction at quality level q at accuracy a. A general observation is that 
with increasing sampling rate u or increasing observation time T, the quality increases, 
due to more accurate information that is obtained about the system's states. In general, 
however, it is not only the total number T x uj of restricting equations (per node and 
per dimension) that controls the quality. For instance, at a given observation time, high 
quality close to Qo.qs = 1 is reached even at lower sampling frequency if the collective 
dynamics is more irregular, cf.. Figure [3^, b vs. Figure [3]3,d. We ascribe this to the lower 
degree of correlation among observations at different times that is required for accurate 
reconstruction numerics in ([g]), e.g. for irregular dynamics, sample points more distant 
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Figure 3. Quality of reconstruction (Q0.95) as function of sampling rate u and 
observation time T. Good quality reconstruction can be achieved even at low sampling 
rates. Underlying system dynamics used for reconstruction in panels (a)-(d) are shown 
in Fig. [1] (a)-(d). Lengths of observations are T = 1 (A), T = 5 (+), T = 10 (0), and 
T = 20 (O)- Each point is the result of averaging over 50 networks. 



in time provide more relevant information increase about the system as they are less 
correlated than closer- by points. 

How does the minimum required observation time scale with system size? Figure [4] 
shows ro.98,0.95,200 5 the minimum length of time of observation required to have at least 
q = 98% of the links accurate in strength to an accuracy of at least a = 0.95, sampled at 
a rate of uj = 200, as a function of A^. The numerics suggest that, at fixed cj, Tq^^.u; scales 
sublinearly with network size A^ for reasonably small < 1 — a <C 1 and < 1 — g <C 1 
(reasonably large a and q). This implies that the cost of observation does not grow 
prohibitively quickly, and that even large networks can be reconstructed by a single 
observation of the collective dynamics. 



3.4' Robustness: Substantial noise and unknown parameters 

Is reconstruction still feasible in the presence of substantial noise? Is it possible 
to find unknown parameters of the local unit systems? In the preceding examples, 
none of the unknown parameters appeared in the dimension of coupling, making the 
problem of reconstructing network topology distinct from that of inferring dynamical 
parameters. In the following example we illustrate reconstructing networks with intrinsic 
unit dynamics that are governed by arbitrary functions with K unknown parameters, 
and where the dynamics are influenced by substantial additive noise ^l . Assuming 
that all Q have a finite variance, we note that an observation of the dynamics of the 
system yields a system of equations linear in the K+N unknowns, which we can solve 
as before. A simple example is a network of Lorenz oscillators [9], where the dynamics 
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Figure 4. Sublinear scaling of reconstruction time with system size, (a) Minimum 
required time of observation for reconstruction (g = 0.98, a = 0.95,cj = 200) grows 
sublinearly (presumably algebraically) with system size N for networks with various 
fixed indegrees K = 10(»), K = 50 (v), K = 100(a). The fit suggests that 
Tq^a,uj ^ N^ ^ where the exponent of scaling 7 ^ 0.5. (b) An algebraic scaling 
(black) fits the data (ir=100, A) best. Data plotted on a bilogarithmic scale. Linear 
best-fits (green dashes) overestimate and logarithmic fits (red dots) underestimate 
reconstruction time. All fits to data from 7V=100 to 7V=2000. 



of each oscillator is given by 



N 



Xi 



CTi {Vi - Xi) + J^ Jij {xj - Xi) + rii^ 



(x) 



J=l 



Vi = Xi {pi - Zi) -yi + rj^l^^ 
Zi = XiVi - /3iZi + r]^^^^ 



(11) 



where the parameters cr^,Pi and /3i are unknown, and chosen randomly from intervals 
where the Lorenz system is known to be chaotic: a^ G [9, 11] ,Pi G [20,35] ^ /3i G [2,3]. 
The noise terms ^l , rf G {x, ^, z}, i G {1, . . . , A^} are independent identically distributed 
random variables (on the discrete simulation time scale of At = 0.001) drawn from 
the standard normal distribution. The network topology and all parameters of the 
system can be reconstructed and the reconstruction method works despite substantial 
interference from noise. Figure [5] shows a successful reconstruction of the network and 
of all dynamical parameters for a network of heterogeneous Lorenz oscillators, where 
the noise amplitude ry = 0.5 is chosen such that it drastically alters the dynamics 
from its deterministic counterpart (black vs. blue curve in Figure 5a). This illustrates 
by example that the theory is insensitive to additive noise and capable of successful 
reconstruction, as desired for generic real-world systems. 
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Figure 5. Reconstructing a network and unknown parameters for a system in the 
presence of substantial external noise, (a) The dynamics of a unit in a network 
of 32 Lorenz oscillators in the noise- free (blue) and noise-driven (black) regimes. 
The network was a realisation from an ensemble of networks with edge connection 
probability p = 0.5. Starting from the same initial condition, the noise-driven 
trajectory quickly deviates due to the chaotic nature of the system. Reconstruction 
of the network topology (d) and parameters (f) with corresponding absolute errors 
(e), (g). (b) shows the actual network, and (c) shows the false positives in the 
reconstruction. There were no false negatives. 



4. Conclusion 

In summary, we have introduced a simple robust method to infer connection topology 
from observations of deterministic and noisy network dynamical systems, where the 
functional form of the evolution equations is known. 

Our method is unique in the following ways. A simple sufficiently long observation of 
the system dynamics suffices to reconstruct the network topology and coupling strengths. 
In many cases, experimental access to the system, say, to introduce 'control signals', as 
in [U [6] may not be possible, and the method proposed here does not require any form 
of system intervention. Furthermore, the type of collective dynamics is not restricted 
to, e.g., fixed points, periodic orbits, synchronous states, or any other specific type 
of motion, cf. O H]. Using entire dynamical trajectories, including transients, for 
the reconstruction, we demonstrate robust reconstruction for a wide range of observed 
dynamical states, from asynchronous chaotic states to transient states towards global 
synchrony. Moreover, we treated the system as a 'gray box', where we know some general 
principles of the system (like the coupling functions) but lack the details, like network 
structure or intrinsic parameters. Given (e.g. experimental) dynamical observation 
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data, our theory provides an explicit analytic solution to the inverse problem of finding 
the network structure. This solution is a direct restatement of the differential equations 
governing the dynamics of the system, and is thus conceptually the simplest possible. 
This simplicity suggests highest attainable quality for such inverse problems. 

The method scales sublinearly with network size, seems robust against substantial 
addition of noise, and thus provides a promising complement to existing reconstruction 
methods. Thus, our method offers a conceptual simplification over other methods that 
make the same assumptions we do, but rely on more complex techniques like copy- 
synchronisation (called auto-synchronization in [6]) or the use of topology estimating 
clone models or control signals [5j. Further, the method suggested here is capable of 
reconstructing network structure from a simple observation of the system's dynamics, 
without resorting to any external intervention to drive the system into or from some 
canonical state, as in [4J. 

Efforts to understand the general interplay of network structure and dynamics 
have yielded several promising approaches, mainly applicable to smaller systems. 
Notable among the forward methods, i.e., methods that predict dynamical features 
from knowledge of network topology are those that study the propagation of a harmonic 
perturbation through a network of coupled phase oscillators [IQl IH] , and methods to 
predict disordered dynamics from structures of strongly connected components in the 
network [12]. Cimponeriu et al. [13j introduce two methods to estimate the interaction 
delay in weak coupling between two self-sustained oscillators from observed dynamical 
time series. Arenas et al. p!4] show that times of synchronisation can reveal the 
hierarchical structure of a network, revealing a connection between synchronisation 
dynamics and topological clustering. In an alternative approach, Memmesheimer and 
Timme p!5l IT6] present an analytical method to design networks of spiking neurons 
that display a required spike pattern. Other inverse methods have relied on stochastic 
optimisation [T7] to fit a model of a network of spiking neurons to an observation of a 
real network to infer its topological parameters. 

Several avenues for further research present themselves. As suggested by previous 
work [4^ i2j, minimising the ii norm, instead of the £2 norm as used here may result 
in more efficient reconstruction of sparse networks [18]. In addition, our preliminary 
studies suggest that this highly overdetermined inverse problem can be reduced to 
an exactly-determined problem by selectively choosing points on the time series to 
ensure that the resultant system of equations is maximally linearly independent. This 
reduction significantly reduces the cost of computation to reconstruct large networks. 
Furthermore, it is straightforward to extend this method to coupled map networks and 
to systems with delay. Finally, recent studies show that a method analogous to the one 
presented here for smoothly coupled systems is capable of reconstructing networks of 
pulse-coupled systems such as integrate-and-fire neurons [19j. 
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